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ABSTRACT 

With  the  advent  of  the  high  speed  digital  computer,  many 
problems  heretofore  considered  unsolvable  for  all  practical 
purposes  are  now  well  within  the  reach  of  the  applied  math- 
ematician.  One  such  problem  is  the  routing  of  a  shin  throuqh 
a  time  dependent  ocean  wave  field,  from  one  point  on  the 
earth's  surface  to  another,  so  as  to  minimize  a  cost  func- 
tion of  the  form  g(x,y,t,u). 

This  paper  considers  a  numerical  solution  to  the  above 
problem.   The  technique  to  be  employed  is  known  as  the  method 
of  steepest  ascent  and  is  attributed  to  Arthur  E.  Bryson  and 
Walter  F.  Denham  [1].   Although  the  computer  program  as  given 
in  the  Appendix  is  written  specifically  for  a  VC2AP3  class 
vessel  operating  in  a  described  area  of  the  North  Pacific 
Ocean,  it  can  be  readily  modified  to  accommodate  any  type 
vessel  operating  in  the  Northern  Hemisphere. 
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I.    Statement  of  the  Problem 

The  basic  problem  to  be  developed  and  solved  in  this 
paper  is  the  following:   A  VC2AP3  vessel,  located  at  latitude 
40°N  longitude  154°E,  is  to  transit  a  time  dependent  ocean 
wave  height  and  wave  direction  field  to  latitude  38°N  longi- 
tude 123°W  so  as  to  minimize  a  cost  function  of  the  form 
g(x,y ,u,t) . 

A  south-polar  stereographic  projection  of  the  Northern 
Hemisphere  upon  a  plane  passing  through  the  circle  of  60°N 
latitude  is  used  to  establish  a  rectangular  coordinate  sys- 
tem, with  the  OX  and  OY  axes  parallel  to  the  projections  of 
the  meridians  of  10°E  and  100°E  longitude,  respectively. 
Then  a  63  x  63  grid  is  constructed  along  the  OX  and  OY  axes, 
with  x  =  y  =  31  defining  the  projection  of  the  North  Pole. 
At  each  point  of  this  grid,  actual  wave  heights  and  wave 
directions  taken  from  the  files  of  the  U.S.  Navy  Fleet  Numer- 
ical Weather  Central,  Monterey,  California,  are  recorded. 
The  map  scale  factor  MSF,  defined  as  the  ratio  of  a  differen- 
tial distance  in  the  OXY  plane  to  the  corresponding  differen- 
tial distance  on  the  earth's  surface,  is 

MSF  =  [973.75+(x-31)2+(y-31)2]/1043.6. 
Germane  to  the  solution  of  the  above  problem  is  the 
elliptical  polar  velocity  diagram  of  Figure  1,  giving  the 
ship's  speed  v  as  a  function  of  the  angle  0  between  the 
ship's  heading  and  the  wave  direction.   R.  W.  James  [6]  gives 
empirical  curves  for  the  speed  of  a  VC2AP3  class  vessel  in 
head  waves  v,  ,  beam  waves  v,  ,  and  following  waves  vf  as  a 


function  of  wave  height.  Bleick  [2]  showed  through  a  least 
squares  analysis  that  these  three  speeds  can  be  represented 
closely  by  the  4-parameter  formula 

v  =  C4(H+C1)  *    exp(-C3H)  , 

where  H  is  the  wave  height. 

Thus,  the  projected  speed  of  the  ship  at  any  point  on 
the  grid  is 

V(x,y,0,t)  =  v(H,0)MSF  , 
where  v(H,0)  is  the  ship's  speed  of  the  polar  velocity  dia- 
gram of  Figure  1,  and  H(x,y,t)  is  the  interoolated  value  of 
the  Fleet  Numerical  Weather  Central  wave  height  data.   If 
we  define  K(x,y,t)  to  be  the  stereographic  grid  wave  direc- 
tion measured  counterclockwise  from  the  OX  axis,  and  u(t) 
to  be  the  ship's  heading  as  measured  counterclockwise  from 
the  OX  axis,  then 

[v  sin(u-K)/b]2  +  {[c+v  cos (u-K) ]/a}2  =  1   or 

V(H,0)  =  b(a  "c  }/[—  cos0  +  /b2cos20+(a2-c2)sin20] 
a       a 

where  0  =  u-K,  a  and  b  are  the  semi-principal  axes  of  Figure 

1,  and  c  is  the  distance  to  the  eccentric  pole  o.   We  note 

at  this  point  that  an  elliptical  polar  velocity  diagram, 

such  as  Figure  1,  must  be  specified  for  each  wave  height, 

and  that  a,  b,  c  are  functions  of  H(x,v,t).   Thus,  we  can 

write  the  equations  of  motion  of  the  ship's  projection  in 

the  OXY  plane  as : 

x(t)  =  Vcos  u(t) 

y(t)  =  Vsin  u(t) 
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where  the  (•)  over  a  variable  indicates  its  time  derivative, 
and  V  =  |v|  =  V(x,y,u,t).   Hence,  our  problem  is  to  choose 
the  control  angle,  u(t) ,  at  each  point  of  the  transit  so  as 
to  minimize  the  cost  function  g. 


II.   Abstract  of  the  Theory 

Consider  a  system  of  differential  equations  of  the  form 

*  i     i  1     i 

x  =  f   (i=l,...n),  where  f  =  f  (X,U,t)  are  functions  of 

2 
class  C  ,  X  is  an  n-dimensional  vector,  U  is  an  m-dimensional 

vector,  and  t,  the  independent  variable,  is  a  scalar.   We 
desire  the  solution  to  the  above  system  of  equations  which 
is  optimal  in  the  sense  of  maximizing  (minimizing)  the  ter- 
minal value  of  some  cost  function  g(X,t) _,  while  simulta- 
neously satisfying  terminal  constraints, 

hi(X,t)t=T  =  0  (i=l,. . .N;  0  £  N  <  n) . 

The  solution  of  a  problem  of  the  above  form  by  the  method  of 
steepest  ascent  is  concerned  with  the  development  of  the  m 
control  variables,  u, (t),...u  (t) .   Essentially  the  oroblem 
is  solved  in  two  parts.   In  part  I  we  are  interested  in 
attaining  admissibility,  that  is  satisfying  the  terminal 
constraints 

hi(X,t)t=T  =  0  (i=l,. ..N;  0  <  N  <  n) . 

In  part  II  we  are  concerned  with  the  problem  of  maximizing 
(minimizing)  the  terminal  value  of  the  cost  function  g.   The 
procedure  is  as  follows:   We  initially  guess  nominal  values 
of  the  control,  U;  in  general,  the  resulting  solution  will 
neither  be  admissible  nor  an  extremal.   Changes  in  the  con- 
trol are  then  generated  which  tend  to  drive  us  to  admissibil- 
ity, that  is,  U(t)  is  replaced  by  U(t)  +  6U(t).   This  nro- 
cedure  is  continued  until  we  attain  admissibility.   The  ad- 
missible solution,  in  general,  will  not  be  an  extremal;  so 


we  once  again  generate  changes  in  the  control  which  tend  to 
drive  us  toward  an  extremal  without  affecting  admissibility. 

The  author  realizes  the  above  is  rather  abbreviated; 
however,  this  discussion  is  intended  to  give  a  brief  idea  of 
what  is  to  follow. 
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III.   Development  of  the  Theory 

Notation 
Let  us  start  by  considering  the  following  system  of 
ordinary  nonlinear  differential  equations: 

x   =x=  Vcos  u=  f  (x  ,x  ,u,t) 

•  2  2   12 
(3.1)   x   =  y  =  V  sin  u  =  f  (x  ,x  ,u,t) 

•  3  3   12 

x   =  z  =  g(x,y,t)  =  f  (x  ,x  ,u,t) 

where  V  =  V(x,y,u,t),  u  =  u(t),  and  g  is  the  cost  function 
to  be  minimized.   We  can  write  this  as 

x1  =  f1(x1(t) ,  x2(t) ,  u(t) ,  t) ,  (i  =  1,2,3) ,  0  <  t  <  T, 
where  t,  the  independent  variable,  is  a  scalar.   It  will  be 
convenient  in  the  development  to  follow  to  represent  these 
sets  of  variables  by  matrices  or  by  vectors,  so  let  us 
define : 


X  = 


xX(t) 
x2(t) 
x3(t) 


;  F  = 


f  (x  ,x  ,u,t) 

2   12 
fz(x"\xz,u,t) 

*3/  1   2    .» 

f  (x  ,x  ,u,t) 


;  U  =  [u(t)] 


where  U  is  a  single-rowed  column  matrix  since  there  is  but 
one  control  variable. 

Equations  (3.1)  then  become 
(3.2)  X  =  F 

The  vector  X  is  called  the  state  vector,  or  state ,  and  the 


one  dimensional  vector  U,  the  control  vector,  or  control. 
We  will  assume  the  control  to  be  a  continuous  function  of 
time . 
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Now  the  problem  as  posed  in  section  I  can  be  interpreted 
as  follows:   We  desire  to  find  a  curve  C,  a  solution  to  (3.2) 
for  some  choice  U,  which  is  optimal  in  the  sense  of  minimizing 
the  terminal  value  of  the  cost  function  g  while  simultaneouslv 
satisfying  the  terminal  constraints 

x1(T)  =  x*,  x2(T)  =  x2f      , 

1       2 
where  xf  and  x,.  are  the  final  values  of  the  state  corresnond- 

ing  to  the  rectangular  coordinates  of  the  terminal  point  in 
the  OXY  plane. 

The  Adjoint  System  of  Equations 
Of  primary  interest  is  the  relationship  between  varia- 
tions, 6u(t) ,  in  the  control  at  any  time,  t  (0    t   T) ,  and 
the  resulting  terminal  variations  of  the  state  variable, 
SX(T).   If  we  consider  (3.2)  in  comoonent  form,  the  varia- 
tional relationships  for  the  state  variables  can  be  expressed 

as  : 

,•1    df1   x    1        Sf1  ,  2    af1  x  3    Sf1  x 

OX   =  =-  OX   +  7T    OX   +  rr    ox   +  ou 

3x         3x         3x         3u 

2  2          2          2 

, -,  0x  j,  •  2    3f  r    1  3f  r  2    3f   r  3    3  f  , 

(3.3)   6x   =  — y   ^x  +  — o"  ^x   +  — T  ^x   +  ^u 

3x  3x         3x         3u 

-•3    3f3  £    1    3f3  ,  2    3f3  .  3    3f3  , 
6x  =  — y   6x  +  — =■  ox  +  — =-  ox  +  ou 

3X1       3x        3x        3u 
Writing  this  in  matrix  form  we  have: 

6X  =  G6X  +  H6U 
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where 


G  = 


8f    9f" 


9f 


ax- 


ax   ax 


3 


3f*   3f 


3x    3x' 


3fJ   3f 


3x    3x' 


2  ' 

3f 
7  i 

8x- 


3x" 


H  = 


"af1! 

3u 
3f2 

3u 
3f3 

3u 

6X  = 


6x- 


5x- 


6x' 


and  6U  =  [6u]  . 

Note  in  the  special  case  where  the  control  U  consists  of  a 
single  variable,  u,  6U  can  be  considered  a  single-rowed 
column  vector. 

If  we  now  introduce  a  new  set  of  variables 

X(t) 

p(t) 


A(t)  = 


ait) 


Lagrange  multipliers,  which  are  undefined  functions  of  t,  and 

T         T 
multiply  (3.4)  through  by  A  ,  where  A   is  the  transoose  of 

A,  we  obtain  the  scalar  relationship 

rn  i       rp    _>.       m    _^ 

(3.5)  A  6X  =  A  G6X  +  A  H6U   . 

Let  us  now  integrate  (3.5)  over  the  interval  (0,T)  to  obtain 


•T  AT0;t 


■T  .T„„:fc 


•T  ,T„.A 


(3.6)  jQ    A    SX   dt  =  jQ    A  G6X  dt  +  jQ    A  H6U  dt  . 
Integrating  the  integral  on  the  left  by  parts  yields 

T  ->  i  T     rT*T->        />TT->-        ^TT-v 

A  6X | J  -  /J  A  6X  dt  =  /J  A  G6X  dt  +  jQ   A  H6U  dt 
Rewriting  this  as 

A  6X\       -    jl     (A   +  A  G)6X  dt  =  j n    A  H5U  dt 
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and  choosing  A  as  a  solution  to  (AT  +  ATG) 6X  =  0  defines  the 
system  of  equations  adjoint  to  the  variational  equations  (3.4), 
That  is,  the  adjoint  system  of  equations  is 

•m         m 

(3.7)  A   =  -AG 

With  this  choice  of  A ,  (3.6)  reduces  to 

(3.8)  AT6X|T  =  ATSX|0  +  /J  ATH6U  dt 

and,  since  we  are  assuming  6X(0)  =  0 ,  we  have  a  relation  be- 
tween the  terminal  variations  of  the  state  variable,  6X(T), 
and  the  variations  in  the  control,  6U.   Now  if  A  is  chosen 
as  the  specific  solution  to  the  adjoint  system  of  equations, 
(3.7),  such  that  at  time  t  =  T,  A(T)  =  (1,0,0),  then  (3.8) 
becomes 

6x1(T)  =  /q  ATH6U  dt 

1  ■* 

giving  the  terminal  variation  of  x   due  to  variations  6U(t). 

J_  1_ 

Similarly,  choosing  A. (T)  such  that  the  i    component  is  1 
and  all  other  components  are  0  yields 

(3.9)  6xx(T)  =  JjJ  ATH6U  dt 

We  may  consider  the  matrix  product  A.H  as  a  vector,  sav  A. (t) , 
when  there  is  more  than  one  control  variable.   Then  (3.9)  can 
be  written  as 

6x1(T)  =  !TQ   Ai(t)  •  6U(t)  dt 

and  for  a  maximum  change  in  x  (T)  we  would  choose  6U(t)  oar- 

allel  to  A.  (t)  ,  for  0  _<  t  <_   T.   Courant  [3]  oaqe  223,  refers 

-*■  i      "*■     "*"   i 

to  the  vector  A. (t)  as  the  gradient  or  x  ,  of  A.  =  Vx  . 

l  *  l 
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Special  Variations 

From  the  previous  discussion  we  have  seen  for  a  maximum 

i  ■*■  a  i 

change  in  x  we  would  like  6U(t)  to  be  parallel  to  vx   for 

0  <  t   T.   This  motivates  the  following  choice  of  a  soecial 

variation  of  the  control , 

(3.10)        6U(t)  =  e1Vx1(t)  +  e2^x2(t)  . 

It  might  be  worth  noting  we  do  not  consider  vx   since  at  this 
juncture  we  are  only  interested  in  satisfying  the  terminal 
constraints 

x  (t)  =  xf,  x  (t)  =  xf  . 


T      ±  1 

Substituting  (3.10)  into  (3.9)  ,  and  recalling  A.H  =  Vx  ,  we 

obtain 

6x1(T)  =  e±    /q  ^x1  •  Vx1  dt  +  e2  /J  V'x1  •  Vx2  dt 
(3.11) 

6x2(t)  =  e±    /q  Vx1  •  Vx2  dt  +  e2  /J  Vx2  •  Vx2  dt 


We  now  define 

0 


Z1^  =  Jj?  Vx1  •  Vx^  dt 


so  that  (3.11)  becomes 

2 
6x1(T)  =  I      e.Z13       (i=l,2)   . 

For  small  perturbations,  the  value  of  T  will  be  changed 
only  a  small  amount  dT,  so  that 

dx1(T)  =  6x1(T)  +  x1(T)  dT   (i=l,2)  . 
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We  can  express  this  relationship  in  matrix  form  as 
r 


z11  z12  k1 


(3.12) 


L 


z21  z22  k2 


~ 

el 

= 

e2 

dT 

_ 

x*  -  xX(T) 
x2f   -  x2(T) 


or  ZE  =  AX.   Let  us  now  digress  for  a  moment,  and  recall  that 
if  we  have  an  overdetermined  system  of  ecruations ,  say  AX  =  B, 
where  A,  the  coefficient  matrix,  is  3  x  2,  X,  the  unknown 
matrix,  is  2  x  1,  and  B  is  a  3  x  1  matrix,  then  a  least 

squares  solution  for  X  is  obtained  by  multiplying  both  sides 

T  T 

of  the  above  equation  on  the  left  by  A  provided  AA   is  not 

singular.   Since  we  have  a  system  of  two  ecruations  in  three 

unknowns,  let  us  apply  the  above  technique  in  reverse. 

c. 


Let   E  =  Z  C  where  C 


so  that  (3.12) 


T  T 

becomes  ZZ  C  =  AX.   Now  ZZ   is  a  2  x  2  Gram  matrix,  so  let 


ZZ   =  A  = 


all   a12 


a21   a22 


Thus  (3.12)  is  representable  as  AC  =  AX,  and  we  are  now  able 

T 
to  solve  for  C.   Use  of  (3.10)  and  E  =  Z  C  gives  us  our  var- 

-> 
lations  in  control  6U(t)  ,  for  0    t  _<  T. 

At  this  point  in  the  development  we  should  be  able  to 

obtain  an  admissible  curve.   We  will  now  consider  the  nrob- 

lem  of  decreasing  the  cost  function  without  affecting  x  (T) 

2 
and  x  (T) .   To  accomplish  this  end  we  consider  a  special 

variation  of  the  form 
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(3.13)  6U(t)  =  a1^x1(t)  +  a2Vx2(t)  +  a3Vx3(t) 
Substituting  (3.13)  into  (3.9)  yields 


6x1(T)  =  I      a. 3! 

j  =  l   3 


i: 


and 


(1=1,2,3) 


dxX(T)  =  6xX(T)  +  xX(T)  dT 
In  matrix  notation  this  becomes 


(3.14) 


X*  -  x1(T) 


x2  -    x2(T) 


\*l   -    x3(T) 


Iz11  z12  z13 


z21  z22  z23 


z31  z32  z33 


xX(T) 

al 

XZ(T) 

a2 

.  3 
xJ(T) 

a3 

_ 

dT 

L 

Since  we  have  an  admissible  path, 

X^  -  XX(T)  a  X2  -  X2(T)  a  0  ; 

3     3 
however,  x   -  x  (T)  is  unknown,  and,  in  general,  is  not  zero 

Once  again  we  are  faced  with  the  problem  of  solving  an  under- 
determined  system  of  equations,  (3.14).   As  before  we  will 

T 
let  a  =  Z  D  where  D  = 


Hd,   so  that  (3.14  becomes  AX  =  ZZ  D. 


|d2 


Now  ZZ   is  a  3  x  3  symmetric  Gram  matrix  so  let 


bll   b12   b13 


ZZ   =  B  = 


b12   b22   b23 


b13   b23   b33 


Thus  (3.14)  is  representable  as 
(3.15)  AX  =  BD 
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Define  d .  =  eg . 0 
1      i3 


(1=1,2,3),  where  e  is  to  be  determined, 
and  3 ■ 3  as  the  cof actor  of  b.3  in  the  determinant  of  B,  that 


is 


dl  =  £[b12  b23  "  b22  b131 


d2  =  e[b12  b13  -  bi;L  b23] 


d3  =  E[bn  b22  -  b12  b12] 


Then,  upon  substitution  into  (3.15),  we  obtain 


x*  -  xX(T) 


(3.16) 


AX  = 


2  2 
xf  -    x    (T) 

3  3 
x^  -  xJ(T) 


0 

=    e 

0 

det  B 

_ 

Thus,  by  choosing  e  negative,  we  can  decrease  x   and  keep  x" 

2 
and  x   fixed  so  long  as  B  is  nonsingular.   That  is 

12  3 

dx    dx    n    ,  dx     ,   _ 

3 —  =  -g —  =  0  and  -g —  =  det  B. 

de     de  de 

T  + 

Now,  having  D  and  a  =  Z  D,  (3.13)  may  be  solved  for  6U  as  a 

multiple  of  e.   In  order  to  limit  the  maximum  change  in  the 

control,  we  define  the  following  norm: 

N  =  {max  |u(t)  |  +  1  Iffil  }   . 
t 

DET  TJMX 

We  now  define  e  = £• ,  where  DELUMX  is  the  maximum 

N 

allowable  change  in  the  control.   Thus  6U(t),  the  solution 
to  (3.13)  is  replaced  by 

DELUMX 


e  6U(t)  =  -- 


N 


5U(t) 
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IV.   The  Numerical  Solution 

The  numerical  solution  of  the  problem  as  posed  in  Sec- 
tion I  is  based  upon  the  solution  of  the  original  system  of 
equations,  the  adjoint  equations,  and  the  variational  equa- 
tions.  For  the  sake  of  simplicity,  a  finite  difference 
approximation  to  the  derivative  with  a  variable  time  step  of 
integration  is  used  to  integrate  the  above  equations.   The 
variable  time  step  of  integration  is  obtained  by  dividing 
the  total  time  as  determined  after  each  path  into  two  hun- 
dred equal  increments.   Initially,  the  value  for  T,  the 
terminal  time,  is  determined  by  the  geodesic  track,  but  this 
is  not  a  requirement.   Any  reasonable  guess  for  T  will  suf- 
fice.  Justification  for  this  method  of  integration  is  based 
upon  a  comparison  of  the  geodesic  route  as  obtained,  using 
the  above  mentioned  technique  and  a  fourth  order  Runge-Kutta 
integration  method. 

Due  to  the  core  storage  limitation  imnosed  by  the  Fleet 
Numerical  Weather  Central,  165,000  octal  words,  a  theatre  of 
operation  is  established  which  is  a  subfield  of  the  63  x  63 
grid  mentioned  in  Section  I.  Having  established  this  sub- 
field  the  wave  field  data  is  then  read  into  core  memory.  We 
are  now  in  a  position  to  proceed  with  the  numerical  solution 
of  the  problem. 

Although  not  a  requirement,  the  geodesic  track  from  the 
point  of  departure  to  the  point  of  destination  is  calculated 
for  comparison  purposes.   As  stated  in  Section  II,  a  nominal 
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control  program  is  initially  required.   This  is  taken  as  the 
angle  of  inclination  of  the  straight  line  connecting  the 
track  end  points. 

Recall  that  in  the  development  of  the  theory  we  chose 
as  particular  solutions  to  the  adjoint  equation,  A,,  A2#  and 
A.,,  where  A.  (T)  has  its  i    component  equal  to  1  and  all  other 
components  equal  to  0.   To  satisfy  this  requirement  we  define 


A1(0)  = 


,  A3(0)  = 


,  A2(0)  =  I  1 

L° 

and  integrate  the  adjoint  system  of  equations  forward  to 
obtain  A,(T),  A2 (T)  ,  A3(T).   The  resulting  matrix 


X1(T)   X2(T)   X3(T) 


y1(T)   u2(T)   y3(T) 


a1(T)      a2(T)   a3(T) 


is  then  inverted  to  obtain  new  initial  values  A,(0),  A2(0), 
A3(0)  so  that  upon  integration  forward  a  second  time  one 
obtains 


0 


0 


A1(T)  =   0  ,  A2(T)  =  •  1  j  ,  and  A3(T)  =  I  0 
0  j  0  [l 

Having  obtained  the  initial  values  of  the  adjoint  vectors 

A,,  A2,  and  A3  such  that  at  the  terminal  time  T  they  take  on 

prescribed  values,  we  will  run  a  path,  check  for  admissibilitv, 

and  if  not  admissible  generate  corrections  to  the  control. 

Simultaneous  integration  of  the  original  system  of  equations 

and  the  adjoint  equations  gives  a  path.   At  each  point  along 

this  path  the  gradients  in  the  x   direction  (i=l,2,3)  are 


20 


calculated  and  stored  for  later  recall  in  determining  the 
changes  in  the  control.   In  addition,  the  elements  of  the  Z 
matrix  of  Section  III  are  calculated.   A  check  is  now  made 
to  see  if  we  have  attained  admissibility.   If  not,  we  oro- 
ceed  to  calculate  the  changes  in  control  which  will  tend  to 
satisfy  the  admissibility  constraint.   Assuming  we  do  not 
have  an  admissible  path,  the  A,  C,  and  E  matrices  of  Section 
III  are  calculated.   Thus  we  obtain  new  values  for  the  con- 
trol and  the  terminal  time  T,  that  is,  U(t)  is  replaced  by 
U(t)  +  6U(t)  and  T  is  replaced  by  T  +  dT .   The  above  proce- 
dure is  continued  until  the  admissibility  constraint  is  met. 

Having  admissibility  we  consider  the  problem  of  de- 
creasing the  cost  function 

fT   Y  H,, 
g  =  J0  eT   dt  , 

where  y    is  a  scalar,  without  affecting  admissibility.   A 
check  is  made  to  see  if  the  terminal  value  of  the  cost  func- 
tion as  determined  by  the  last  admissible  path  is  less  than 
that  determined  by  the  previous  admissible  path.   If  it  is, 
we  proceed  to  calculate  the  changes  in  control  as  defined  in 
Section  III,  so  as  to  decrease  the  cost  function  g.   If  not, 
we  decrease  the  maximum  allowable  change  in  the  control  at 
any  time  t  by  one-half.   A  check  is  then  made  to  determine 
if  the  maximum  allowable  change  in  control  is  less  than  0.01. 
If  it  is ,  a  solution  has  been  obtained.   If  not,  proceed  to 
calculate  the  changes  in  control  so  as  to  decrease  the  cost 
function.   Assuming  we  have  no  solution  we  calculate  the  B, 
D,  and  a  matrices  of  Section  III.   This  enables  one  to  obtain 
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new  values  of  the  control,  that  is,  U(t)  is  replaced  by 
U(t)  +  6U(t).   A  new  path  is  run  which,  in  general,  is  not 
admissible;  so  we  again  drive  to  admissibility.   This  proce- 
dure is  continued  until  a  solution  is  obtained. 

It  should  be  noted  the  values  used  in  the  program  for 
the  maximum  allowable  change  in  the  control  and  the  stopping 
criterion  are  quite  arbitrary. 
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V.    Conclusions 

The  value  of  this  study  lies  not  in  the  solution  of  the 
particular  problem  considered,  but  in  the  fact  that  given 
any  realistic  cost  function  which  is  a  function  of  the 
stereographic  grid  coordinates  x  and  y,  the  control  U,  and 
time,  one  can  obtain  an  optimum  solution  with  only  slight 
modifications.   Such  a  cost  function  for  a  cargo  type  vessel 
might  reflect  costs  due  to  storm  damage,  failure  to  meet 
scheduled  commitment  dates,  and  spoilage  of  cargo  due  to  ex- 
tended time  at  sea,  to  mention  but  a  few. 

On  the  basis  of  the  results  obtained  the  following 
observations  seem  justified: 

(a)  It  was  shown  in  Section  III  that  a  neces- 
sary condition  for  the  solution  curve  C  to  be  an  ex- 
tremal is  that  the  determinant  of  the  B  matrix  at 

t  =  T  approach  zero  as  track  iteration  continues. 
In  the  case  where  we  chose  y  =  0,  that  is,  minimum 
time,  we  saw  that  this,  in  fact,  was  the  case.   In 
the  second  problem,  where  we  chose  y  =  0.01,  the 
determinant  of  the  B  matrix  was  reduced  from  a  value 
on  the  order  of  5,000  to  0.9.   The  conclusion  to  be 
drawn  here  is  that,  although  the  solution  curve  C 
is  not  an  extremal,  it  very  closely  approximates  one. 

(b)  Defining  the  initial  control  as  the  con- 
stant function  u(t)  =  V  (0  _<  t  _<  T)  ,  where  V  is  the 
initial  angle  of  departure,  leads  to  convergence  to 
within  fifteen  miles  of  the  point  of  destination. 
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Thus  there  seems  to  be  no  problem  as  to  how  to 
initially  choose  the  control  so  as  to  obtain  con- 
vergence. 

(c)  The  simplified  integration  routine  used 
is  feasible  and  yields  satisfactory  results. 

(d)  The  requirement  of  storing  an  entire  con- 
trol program  and  a  complete  set  of  functions 

T 
A.H   i  =  1,2,3  for  any  one  iteration  is  a  marked 

disadvantage.   Due  to  the  large  storage  recruire- 

ments  imposed  by  the  establishment  of  the  wave 

fields  a  relatively  small  amount  of  core  memory 

remains . 

(e)  An  analysis  of  the  results  obtained  show 
a  substantial  decrease  in  the  wave  heights  for  the 
minimum  time  route,  as  compared  to  the  geodesic 
track,  with  a  decrease  in  the  total  transit  time 
of  approximately  2%  hours.   The  weighted  time 
route,  on  the  other  hand,  shows  only  slight  de- 
creases in  the  wave  heights  encountered,  as  com- 
pared to  the  minimal  time  track.   This  leads  one 
to  conclude  that  if  the  objective  is  to  limit  the 
wave  heights  encountered,  a  function  other  than 
exp  yH  should  be  used.   This  is  more  readily  seen 
if  we  expand  e   ,  that  is, 

eYH  ,  1    +    yH  +  (YH)i+  <pl+    ... 

For  small  values  of  y  and  nominal  wave  heights,  H, 

the  dominant  factor  in  minimizing  z  =  JQ  e    dt  is 
time.  ~  . 
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APPENDIX 


COMPUTER  PROGRAM 
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pijOGUAM    TFSl      (  I\iPiir»0UT^'iT»TAPE4) 

i   ruFNSION    '^TG(12)  ill  (200)  »DELI  1(200 )  fGPAOX( 200)  , GRADY  (200) , 
1     rPS(20l) tHT(201) »A7(20l) tAL(20l) »WA(20l) ?WE(39H9) tWn(63t63> 
m<I.AM(3)  »TI  AM  (3)  tXMU(3)  ?XSIGma (3)  ♦GRaD7(20  0)  #TMH(2) 

COMMON  vf Yt A,R,CC«H,CK,SKf KXtKYtKT.FQXtFOYtKCON 

COMMON /I  l/XhT(8443)  *C^* (8448)  tSNK  (844  8)  »T 

rOMNON/L2/HXtHY 

C-iMMON/LS/CAPV,  CAPVX  ,rAPVY,CA»»\/U»ui 

j  OMM0N/Lf»/XK  »XLG«  XLT 


t  (.Mil  t/ALFNCE     (*E  ( 1  )  two  (  I  f  1  ) 
IMTEPFR    OTG 


) 


FIELD  IM  OUo  THFATFP  OF  INTFOFST 
JMINX=0  %    KTX=0 


*F     I'll  I.    FIRST    EMARl.TSH    THF     <.-AVF 

T  r  A  U  X  s  0     *     0  T  G  X  s  0     *     I  M I  M  X  s  0     5 
pfrr.j)      ^MF/NSTOM^     AMp     MINIMUM     INDICES     OF     ApRAYS 

PF AU     5 • H  V • K Y , K  r  «  T  M I M  , JMI M9 CR  TT 
-.    FORMAT     (     t>T2*F4,0     ) 

t'^IiMT     3tKX»KYtKTtIMIMt  JMTN     fr.PlT 
3    FOpMAT(3^HnpTMFMSI0N    01-     WAVE     FTELl)    ARPA  YS  =  (  I  3  ♦  1  H  «  T3  «  1  H  t  I  3  t  1  H)  4  X5h  T 
lMlMrT3t4xbHjMlfvlaI3»4x5MCRlTsF4.0) 

l'>ys3?-TM'Ikl 

joY=3?-jMl  M 

FOX=IOX 

f OYsjOY 

KCOMsKX»,<Y  +  KX 

y'-F  ANsFl  OAjF  (KX*1  )  /?. 

Y^EAMsF  |OATF(KY*l)/?. 

u  vx  =  XMFAM-2.ni 

POy- yMF  AM„?  ,  !1  1 

oi  fa-  dig  array 

f,()    ^     Js  l  «  k  T 

^   i ,  t  g  ( i )  s  n 
WF  ii.j    octal    PATF.TIME    groups    of    thf    Time    ^fpIfs   mempfrs 

PF  AU     )lt     r»Tf»(l  )  *  ITAU 

J  I    y '">RfaT (o20»HXt 12) 

IF  (OTGX .FO.PTG ( 1  )  .A MO. T  TaUX.FO. IT  All* AND- T M I  NX , FO. IMI M« ANP, 
I  .i^T«Mx.f  O.JMTw.AMp.K  rv  .f-o.KT)  750 #625 
SF T    aWPaY    PTG    to    SELFCTEO    TIme    SFPTE5 
frt>  -,     J   isKT-J 

I    '         /      I  S  1  *  J  J 

CALL  DTPh  (DTG(I) »lTAUfDTG(I*l)  ) 
f   ommttmhF 

pi.r.cy^     THF  UNPACKED  ARRAYS 
PPwlND  4 
'  : )  /on  K s  I  v  K  T 
kkk=(  (  k-1  >*KY-1  )»KX 
>)0  h^fl  KKsl  »2 

MOsy 

K«i  =  KK-l 

r ALL  F IFI  O^X(wE*nTG(K)«KW,NG) 
[F  (\R-]  )  90«f*9990 
M  w  c  w  j  fsi  f  S  ]  ,  0  T  G  (  K  ) 
9)     FORMAT ( 14H0FIELD  OF  PTG=O20 t 1 X31HNOT  FOUND  ON  INPUT  TAPF4,  ARoRT/) 

GO  10  *i 
9n  oo  to  (9;>»<?4)  kk 

9?  DO  H  Is)  ,KX 

I  IsKKK  +  I 
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iif-  l.x  =  !-Tn* 

DO  tt  Jrl.KY 

JJaK.X*J*II 
DFLYsJ-JOY 

ROOT*SQPTF     (    DFLX*DFI_X    ♦    DELY*DELY    ) 
ARR*WD(T*TMlMtJ*J^IM)/5,729577951 
COS*C0SF(ARG) 
STMxSTMF  (  AR(S) 

CSK (JJ)=(-DELX*COS-DFLY*SIN)/POOT 
«  5NK  (  JJ)  =  (OFL  X*STN-DFLY»COS)  /ROOT 
c,a  70  £  SO 

94  i.»n  in  1  =  1  »kx 

TT=  KKK  ♦  I 
D'1  10  J=1,KY 
JJ  =  KX*J  ♦  IT 
in  XHT ( JJ) swn (I*IMIN, J+JMIN) 

650  CONTINUF 

IF(SENSF  SWITCH  l)(S10f700 
MO  print  lino.  KtDTG(K) 
1100  FOPMAt  (lHn»I2f?H.  tO?fl) 

00  t>20  Jsl .KY 
6^0  PRINT  l?00t  (  wO  (  I  +  TMIN»J*JMTN  )#I  =  1»KX  ) 
1200  F OpMAT (  lX»26F5.1/  ) 
700  c^NTTNUF 

PFwlND  4 

DTGX  =  DTGU)  *  ITAUX=TTAU  $  TmtNX  =  TMTm  %    JMINX=JMTN  $  KTX  =  KT 

00  750  T=1,KX 

DO  750  J=1,KY 

JJsKX*J  ♦  I 

RT«I*IMIM  $  Rjsj+JMJN  $  KFPsO 

C-^LL  5F.AI  AMP  (PI,RJ,LAMD,KFP) 

IF  <IVER)646,64H,646 
*4^    PRINT    647       $       GO    TO    B3 

but    FORMAT (26H0RI ,RJ  INOFX  ERPORSt  APORT/) 
64fl  TF  (LAMP) 649,750,649 
644  00  749  K=1,KT 

KKK= ( (K-l ) »KY-1 ) *KX     ♦     JJ 

749  XHT (KKK) =20.0 

750  COM!  IMIjF 

C     *»E  SHOULD  NOW  HAVF  OUR  WAVE  FIELD  OATa  STORED  IN  HFMMOPY 
r     Wfaf)  HOUR  OF  PFPaRTUPF  AND  COORDINATES  OF  TRACK  FMD  POINTS 
hFAD  l3.HR,XLGlfXl.Tl.XLG2tXLT2 
1 3  FO«VAT(F3.n,F6.1fF5.1 ,F6.1 ,F5.1  ) 
PRINT  14,HPtXLGltXLTltXLG2»XLT2 
1*4  F0PMAT(21H  ROUTE  OF  SHIP  BEG  T  MSF4 , 0  ,  31 H  HOURS  AFTFR  TIME  SFPIFS  OR 
1IGIN/19H   FpOM  LONGITUDE  =  F6.1»16h  AMD  LATITUDE  =  F6.1»/19H     TO 
2  LONGITUDE  ■  F6.1tl6H  AMD  LATITUDF  =  F6.1) 
r    wF  'MLL  MOw  CONVERT  COORDINATES  OF  TRACK  END  POINTS  TO  GRID  VALUES 
A»G»{XLGl~10.)/57 ,29577951 
COSLGlsCOSF(ARG) 
SINLGls^INF(ARG) 
ftPGs(X|  O?-10.)/57. 29577951 
QDSLG2*COSF (ARG) 
STMLG2=SINF( ARG) 
MUCScXLTl/5  7.?9b77951 
C05LT1 =COSF ( ARG) 
SINLTlsSlNF < APG) 
ARG=XLT?/57.?9577951 
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♦    nFl   Y»|)FI  Y) 


CnSLT2  =  CnSf-  (APG) 

«;TML  T2  =  MNF  (  A«G) 

k«lsH .pnSfrCOSLTI /(l .+SINLT1 ) 

xi=pp1»coslgi 

Yl=PPl*SlNLGl 

P»2s31  ,2n5»Cost.T2/(l  .♦SIMLT2) 

x?  =  pp?*rnSLC-P 

v  ->  =  pw2*S  j  M|  o? 

v^T»wT=y 1 ♦  Ff'v 

YSTAPT=Y I *FOY 

vF  M")=X2  +  F0X 

r  fi  4-  HTF    GFOOFSTC    kniJTE 

M  sSINUT?»COSLTl*SINLfii|     -    CO<;LT?»«;TN|  Tl  *STNlG? 
tMs    -    STNLT2»C0SLT1*C0SLG1    ♦    CU5LT?*SINLT1*C0SLG2 
t  n=     (SIMLG2*C0SLG1       -    C0SLG2#SINLG1     )  *C0<;|  Tl  »C0SLT2 
POOTsSGRTF  (FL*FL     ♦    FM*FI*    ♦    EN*EN     ) 
H  =EL/R0OT 

FMaEM/POOT 

t  ri  =  EM/pnoT 
HF LX=X?-X 1 
i  H.  Y  =  Y?-Y  1 
S\?  =  SQPTf-  (OFl.X#nFLX 

JF  (XLP^-Xl-^J  )  20,21  ,20 
pn     'A^GsABSF  tfN/*£.4l  ) 

-•."C  =  ASIIjF  (  ARG*S]  2)  /apg 
21     x=xSTART 

r=YSrAPT 

TsHR/24. 

»  TKSrO.O 

CALL    TF  PP 

nT  (  1 )=h 

T  F  S  (  1  )  =  X  T  F  S 

CALL    ANGLF 

ft 7  ( 1)  =xi .G 

Al    (  l)=Xl   T 

-  .A  (  I  )  s  X  K 

MK  =  1 

XX=XSTAPT 
fY=YSTfiPT 

s  s  =  0  .  0 
I.  P  =  0 

STFP=1 ./*. 
00  ^n  K=2f200 

Cr>SPs(F0Y-Y)*FN/31.205  ♦ 
MNP*(X-FOX)*EN/31,205  - 
CALL  TF.RP 
ShlP 

VOF  PT  V 

*   capv«cosp*«;tep 

♦  CAp\/*sInP*STEP 

♦  capv*STEP 
IF (S-APC) 3ht  34»34 

34    PATs(ARC-SS)/(S-SS) 

T=RAT#STFP  ♦  T 
XTFSsRAT*STEB  +  XTFS 
<=(X-XX)»PAT  ♦  XX 


EM 
EL 


CALL 
CALL 

x  =  XX 

ysVV 

s=ss 
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-  MK) 16,16,15 


Y  = ( Y-YY ) owAl  ♦  YY 
GO  TO  IB 
3b  X*  =  X 
YY  =  Y 

ss=s 

T=T  ♦  STFP 

XTFS=X1FS    ♦    STFP 

IF  (  (K-l  )«?_/*  *  1 
]5    NK'sNK    ♦     1 

C*|_L    TEPP 

HT (NK) sH 

TFS(NiK)sXTFS 

C^iL    ANGLE 

AZ(NK)=XLG 

AL (NK ) =XLT 

wa (NK)*XK 

IF(S-APC) 16,41,41 
1<S     \r  (AHSF  (X-XfFAN)  -pox)  37»524t5?4 

3  7     IP  (AHSF  (  Y-Y^'F  AN) -POy)  40*524,  524 

4  0  CONTINUE1 
4]  P9INT  17 

)7  FORMAT(////3oX]4HGEOOESTC  POUTE) 

PPINT  lfl 
lb    FOMMAT(  5X4h0AyS6X5HL0NGI4X5Hl  ATI-5X4HWAVE5X14HWAVE  0IPECTI0M/2X9H 
1  OF  T*?AVF'L4x5H-TUnE4X^HTlJnF5X6HHElRHT6XiOHFROM  NlOPTH/) 

PPT  NT  Jgf(TKS(K),AZ(K)»AL(K)tHT(K)»WA(K)«Kal,MK) 
icj  fowmAj  (  F9,2,FU.ltFa#l,F10,2tF13,Q  ) 
r  we-  WILL  MOW  6UFSS  AN  INITIAL  VALUE  OF  THf  COmTPOl  U(T)  AMD  RUM  A  PATH. 

r)FLTAX  =  XFNn-XSTAPT 

uFLTAy=YFNP-ySTAPT 

r^  (i)F.LTAV  .EO.  0-0)60  TO  530 

ASGsHEl  TAY/OFLTAX 

A|  F A=AT AmF ( APG) 

IF(DELTAx)53l,530f532 

53?  TF  (t)ELT/«Y)  533*534,534 

533  v=8.2*31P  +  ALFA 
GO  7  0  535 

534  v=  ALFA 
GO  TO  535 

531  v  =  3  .  1  4  1  5  P  +  ALFA 

GO    TO    535 
530    IF(UELTAY)536,537,538 

536  \z-=4.7124 
GO    TO    535 

53*    v=1.57o« 
00    TO    535 

537  PPTNT    539 

53P    FOk?MAT(5x*HST0P) 
r-0    TO    P3 

535  CONTINUE 
MWsq 

GAwA=n.n 
541     L««l 

PKX»FL0ATF(KX)-1.01 

«KY=FLOATF (KY) -1 .01 

L«n 

HF|   UMXrn.4 
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r'»i'(HD  =  .ln.(i 

?ni_fjsionoo, 

TAlisXTFS 

nn  boo   t  =  1 • ? o o 

t-'iii    li  (  I  )  =v 

50]     TI.AM  (1  )  =1  .0 

xm(i(  i  jso.n 

X SIGMA  (1  )=(..() 

1I.AM(?)  =0.0 

XMU(2) =1 .0 

XSTb-lA  (?)  =11.0 

1 1.  A  M  (  3 )  =  0  •  0 

KMll(  3)  sO.O 

XSIGMA  (3) =1  .0 

X=XSTAPT 

y=YSTAPT 

$TF.P=TAU/?O0. 

T=HR/24, 

[•0    50?     Tsl  ,20  0 

CALL    TF.PP 

aPG=GAMA»M 

CALL    SHIP 
w    =    1 1  (  I  ) 
CM.L    VDFRIV 
no    550    J=l,3 

XI   AM  (J)  sTL  At/  (J)-(CAPVX*COSF  (II  (I)  )*TLAM(J)  ♦CftPVX*SlNF  (U(I  )  )*XMII(J) 
1     ♦  GAMA<*F*PF  (  APG)  »Hx*xSTGMA  (  J)      )*STFP 

xmu(j)«xmu(j)-(cap\/y*cosf  (U(t)  )*tlaM(  j)  *capvy»sinf  (Ud)  )»xmij(j) 

1     *GAMA*FXPF-  (ARG)  #HY»X5J.GMA  (  J)      )  *STFP 
55  o    Tl  AM(  J)=XLAM(  J) 

XsX*CAPV*CoSP (U(T) )*CTFP 

Ysy*CAPV#STNF  Wit  I )  ) ttSTFP 

IF (X     ,L  T,     ?.0l)GO     TO    5?4 

IF  (  Y     ,LT.     ?«Ol)Gn    TO    5?4 

ir  (X     .GT.     PKX)GO    TO    5?4 

IF(Y     .GT.     PKY)GO    TO    524 

T=T+STFP 
50?    CONTINUE 

T'V|)(  1  )bXMIJ(1) 

T!HI(?)  sXMU(2) 

onjDFTs    XLAM(l)  *XMIJ(2)-XLAM(2)  *XMii(l  J 

1LAM( ] )sxMu(2)/AD iOFT 

Xmik  i) = - x M t j  ( i ) /ADJDET 

x^7GMA( 1 ) sO.O 

n   AM(2)s-XLAM  (?)/A|)JOFT 

ximi  (?)  -XI  AM(  ]  )  /AOJOFT 

X^IG^A (2) sO.O 

TLAM(3)  s(XlAM(2)»XMU(3)-xLAM(3)*TM|j(?)  ) /ADJDFT 

XMH(3)s(XLAM(3)»TMU(  1)-XLAM(D*XMU(3)  ) /ADJDET 

XSI6MA (3) =1 ,n 

X=XSTAPT 

VsYSTAFT 

7  =  0.0 

T=HR/?4. 

XTFSsq.O 

211*0.0 

7l2sO,0 

713=0,0 
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5M 


503 


51  1 


722*0. 0 
Z23s0#0 

733=0.0 

DO  503  1=1 ,^00 

CALL  TFPP 

«Pf,sfiAMA#H 

HT(1)=H 
TFM  I  )=XTFS 

CALL  ANGLE 

«7(  I  )  =XI  G 

-U.  (I)=X|  T 

in  A  (  I  )  *  X  K 

CALL  SHIP 

W  =  I »  (  I  ) 

CALL  VDFPI^ 

G»ADx(  I)sTLAM(i  )*  (CAPVU*COSF(||(  I)  )  -CaPv#SIMF  (ii(  I )  )  )  ♦  X  M « j  <  1 )  « 
1  (CAPVU*SIi\'F  OKI)  )  ♦  CaP^*C05F  (>J(  I  )  )  ) 

r,«ADY(]  )  =Ti _AM(p)»(CAPVU*COSF(tKI)  )  -CaPV«MNF  (IJ(I)  )  )  ♦  XM!K2)» 
1  (CAPvi'*5IMF  (U(  T)  )  ♦  CaPv«COSF(U(I)  )  ) 

GPA07(I)=TLAM(3)»(CAPVU*COSF(U(I) )-CAPV*SINF(U(I) )  )  ♦  XMIK3)* 
1  (rAPVl'»SIivF-  di(T)  )  ♦  CAPV*C05F(U(I)  )  ) 

7.11=ZH*GPADx(  I)#OPAny  (I)«5TFP 

7l2sZl2*GPADX(  I)*GRADY(n»STFP 

Z13  =  Z13*oPaL>X  (T  )*GPADZ  (I)»STpP 

7?2  =  z22  +  GPi\DY(  I)»GPADy(  I)#STFP 

/?3=Z23*RRA0Y ( I ) «GRADZ ( I ) »STFP 

733  =  733*r,PA07  (I  )  *o«AD7  (I  )«STFP 

DO  551  J=l,3 

*LaM(  J)sTLa*"  (  J)-(CaPVX»C05F  (ll(  I)  )*TLaM(J)  *CAPV  X*S  TNF  (  U  <  T  )  )»XMH(  J) 
1  +GAMA*FXPF ( ARG) ttHX<*XSIOMA ( J)   )«STFP 

XMU  (J)bXMH  (J)  -(CAPVY«CnSF  ( l.l  (  J  )  )  #TL  AM  (  J )  ♦CAPVY*SINF(U(I>  )*XMH(  J) 
1  +GA^A«FXPF ( APG)*HY*XSIGMA ( J)  ) «STFP 

T I  A  M  (  J  )  =  X  L  A  *  (  J  ) 

x=x  +  CAPV*COSF  (U(  T)  )  **STFP 

Y=Y+CAPV«SINF(U( I) )*5TFP 

?=Z  ♦  pxpF (APG) *STEP 

1=T*STEP 
xTFS=XTF9*STFP 

CONTINUE 
XOOfaCAPV»COSF (IK200) \ 

Yi,)OT  =  CAPV»?lMF(iJ(200)  ) 

APGsGAKAoH 

ZDOTsEXPF ( APrt) 

CALL  Tfc'PP 

^r (201 ) =H 

TFS(?01 )rXTFS 

CALL    ANGI  F 

A7(201)sxLG 

**l    (20l)=XLT 
*.  A  (  2  0  1  )  =  X  n 
XOlFFsXF^n-X 
YiKFF  =  YFMD-Y 

DFLX«X-F0x 

m(tly  =  Y-FOY 

~|SF=(DELX»*2  +  DFLY<»»2*<?73.75)/1043.^3R743 

FPR0P  =  X0IFF»#2    ♦     Yf)IFF##? 

IF     (    FPPOP-MSF»MSF*CPIT»CPIT*0.000023fcf    )     505»505»511 

Al  l=7ll**»2  +  Z12»»?*XOnTo»? 
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^Oh 


Sn7 
510 


512 

544 


54^ 
5ft* 


Al2=ni*71?  +  722*712*XDOT*YUOT 

«  22»  7 1 2**2  *Z22**2  ^00^*2 

RETsAl 1*A22-A12**2 

CI = (A22*XDTFF-A12*YDTFF> /PET 

C?=( All*YDTFF-Al2*XPIFF)/DET 

F l=Zll*Cl*7l2*C2 

F?sZl2»C] *722»C2 

nF|.  TAllsr.l  «VDOT*C2*YDOT 

T  AU  =  7  Ali*nF|.TA(i 

■in   soft   k=1»2oo 

OELU  (K)sFl*GRAPX(K)*F2*GRA[)Y(K) 

1 1  { K ) =U(K)  ♦oELU(K) 

N  =  N*  1 

IF (N    tGT.    25) GO    TO    50  7 

SO    TO    501 

PRINT    51 f) 

F  mwmaT (/SX [4HNO    CONVERGENCE) 

RO    TO    A3 

CONTINUE 

L  =  L.*l 

IF (NNtEO.l )     RO    TO    544 

TF(IA|j     .IT.     TAijOLOjOO    TO    512 

OFLUHXaPELUMx/2. 

IF(DELUMX  .LT,  n.01)  GO  TO  514 

TAUOLOsTAII 

RO  TO  54"S 

IF  (Z.LT.ZOI  [U  RO  TO  545 

DEI  UMX=OFIJIMX/2. 

IF (DFLuMx.LT.n.Ol )     GO    TO    514 

701  0=Z 

♦  Zl2*7l2 

♦  712*722 

♦  Z?3*Zl2 

♦  722*722 

♦  Z23»Z22 

♦  Z23*Z23 


Kll=7ll*7ll 

Hl2=7l?»7ll 

*13SZ13*ZH 

H22=Z1^*Z12 
P23=7l3#7l2 
^33=713*713 


713*713 
723*713 

Z33*?13 

Z23»723 
733*723 
733*733 


xoot*xoot 

Xl)OT*YOOT 
XDoT*ZOOT 
YOOT*YnOT 
Y0nT»Z0OT 
700T*ZOOT 


i)  =  r11*(F?22*b33-h23*h23)-h12*(r12»m33.h23»R13) 
1  ♦ni3*(H]2*b23-P22»Hl3) 
i>1  =  (H13*'i22-H12*F123)  /D 
U2=(Pll*H23-Hl3*Rl2)/0 
L-3-(bl2*bl2-Bll*B22)/L) 
Hl=Zll#Ol  +  712*02  ♦  7l3ttD3 

723*03 

733*03 

♦  700T*03 


Z22*02  ♦ 

Z23*02  ♦ 

YO0T*02 


♦  R2*GRA0Y(D  ♦  B3*GRADZ(I) 


R2=Z12*D1  ♦ 

H3=7l3*Dl  ♦ 

R4=XD0T*0l  + 

RIG= 0.0 

00  515  Isl,200 

DFLU{I)bB1*GRADX(I) 

RR=A8SF  (  DELI  I  ( 7)  ) 

IF(  HB  .GT.  RlR  )BIG=BB 
51?  CONTINUE 

R^  =  OELUMX/(  8TG  ♦1.0*ABSF(B4)/TaU  ) 

DO  520  I=lt200 

)FLU(I)=RR*DELU(I) 
5?o  Utl)sU(i)  ♦  OELU(I) 

IF (L  ,GT.  25)G0  TO  521 

GO  TO  5ni 
5J4  CONTIMUF 
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IF  (  NN  ,GT,  n  ) GO  T^  54  0 
WRIIMT  5?? 
5??  FORMAT  (////2RX1AHMIMTMUM  TIMF  ROUTE) 

PPINT  543»0 
^4  *  FORMAT  (//20X7HDET  R  =F6.?) 

P^TNT  19,  (  TFS(K) »A7(K ) »AL(K) »HT(K) ,WA (K) »Ksl t20] »10  > 

r,  A  M  A  s  0  .  0  1 

MNJsl 

GO  TO  S4  1 
S«+(^  POINT  S^r5 
S«>>  FOPMAT(////?«XiQHWFir,HTEn  TIMF  ROUTE) 

RRINT  543, n 


POINT  1R 
fOlNT  19, 

PO  to  P3 
PRINT  b?3 
FORMAT (5XRHN0 
PO  To  83 
P  R  T  N  T  b  2  5  •  X  *  Y 
F  OPH4T  (RXiffOHSHTP 
1  Sx3hY  splO.S) 
H3  STOP 
F.MO 


R?] 
5?3 

525 


(  TFS(K) »A7(K) «AL(K) tHT(K) »WA(K> tK»l»201*lO  > 


SOL  N) 


TS  MOW  OUTSTOE  THEATER  OF  0PFPATI0N«*X3HX  sFio.Si 
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SUBROUTINE  TFRP 

DIMENSION  C(4) ,P(4) »0(4) ,PX(4) tQY(4) ,HP(4) •CP(4) tSP(4) ,HPX(4) f 
1HPYI4) ,CPX(4) ,CPY(4) ,SPX(4) tSPY(4) ,H0(4),CD(4) ,S0(4),HS(4) ,CS(4), 
2SS<4) ♦HT(4,4),CT(4,4),ST(4,4) 

COMMON    X,Y,AfR»CC«H»CK.SK»KX.KY,KT,FOX,FOY,KCOM 

COMMON/L 1/XHT(A44R) ,CSK (fl448) ,SNK(*448) »T 

COMMON/LP/HX|Hy 

OOMMON/L^/OKXtOKY 

L=XF I XF (T) 

TT=(-INTF(T) ♦T)»2.-l. 

TP1=TT*1. 

TMl=TT-l . 

T?M=TP1#TM1 

IF  (L-KT*3) 1 »1 ,4 

1  IF(L)?.?.3 

2  K4=3 
TM3sTT-3. 

C  (  ] )=TMl#TM3/fl. 
C(2)=-lP]#TN3/4, 
C(3)sT2M/B. 
r,o  in  16 

3  K4s4 

Sa  (3.*TT+2t  >*TT-9, 

F=-4,«TT*0 

C ( 1 )s-T?m#TM1/16. 

C(2)aG#TMl/l6. 

r  (3)s-F*TPl/16, 
C (4)sT2M#TPl/lft. 
GO  TO  IS 

4  IF  (L-KT)6»5»5 

5  K4  =  l 

L=KT-1 
C  ( 1  )  =1  . 

GO    TO    If 
h    Gs(3,*TT*2,)»TT-9. 
K=-4,»TT*G 
CH  >e-T?M*TM]/16, 
TF(L-KT*1)  11,10,5 
in    k<>=2 

C(2)s(G»TMU  (T2M.F)#TP1  )/]6f 
GO    TO    15 

11    r;(?)=G*TMi/ifc. 

IF(L-KT*?)3f 12,10 
13    K4=3 

C(3)s(T2M-F)*TPl/16. 

15    |.  =  L-l 

1^    M=xFT  XF  ( X) -2 

NsXF IXF (Y)-2 

XXs(-lMTF(X) ♦X)»2.-l. 

YY=(-IMTF(Y) ♦Y)#2.-l. 

XP1=XX+1. 

xm1»xx-1 . 

yP]=YY*l, 
r^=YY-l. 
X2M«XPl*XMl 

y?m=ypi»ym1 

P(1)=-XM1«X2M/16. 
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]  7 


1« 


11 


20 


21 


2? 


2.3 


?S 


?* 


P(2>  = 
P{3)  = 

PI4I* 
0(1)  = 
0(2)  = 
u  ( 3 )  ■ 
0(4)s 
PX  (4) 
PX  (1  ) 
ox  (2) 
PX  (31 
OY  (4) 
(jY(1) 
0Y(2) 
'0  Y  (  3  ) 
DO  27 
HP  (K) 
ro(K) 
SP(K) 
H&X  (K 
nPY  (K 

C  P  X  (  K 
C  r>  y  (  K 

SPK  (K 
S  P  Y  ( K 

K*  =  ( 

DO  23 

HO(  J, 

C  0  ( J ) 
SD(J) 
HS(J) 
CS(J) 
SS(J) 
JJaJ* 

no  23 
r  I  =  I  * 

*T(I» 
CT(Ii 

7  T  <  I  . 

no  25 
no  25 

HO(I) 
CD(I) 
SD(I) 
hS(U 
CS(I) 
SMI) 
C^WTI 
no  27 

HP(K) 
CP(K) 
SP(K) 
HOX  (K 
CPX  (K 
SraX  (K 
mPY  (K 
CPY  (K 
SPY(K 


( (3,»XX*2. 
-XX»XX/8.* 
XPl*X2M/16 
-YM1«y2m/1 
(  (3.«YY*2. 
-YY#YY/8.* 
YPl»Y2M/l* 
= (0.375#XX 
=0.5*XX-PX 
=(1 ,12&*XX 

=-o,5*xx-p 

=(0.375*YY 
=  0.S*YY-(jY 
=(1.125»YY 
*-0.5*YY-Q 

K  =  ] «K4 
s  0.0 
=  0.0 
=  0.0 
)=0.0 
)=P  .0 

)  =  0  .  o 

)=0.0 
)  =  0  .  0 
)=0,0 

(K*L ) »KY*N 
J=l#4 

-0.0 

=  0.n 
=  0.0 
=  0,0 
=  0.0 
=  0,0 
KX*KK 

1  =  1  »4 
J  J 

J)=xpt  (in 

J)sCSK(II) 
J)sSNK(TI) 

1  =  1  #4 

J=]  »4 
=Q(J)#HT(I 
*0( J)*CT(I 
■Q(J)*ST(X 
=P( J)*HT(J 
=P(J)#CT(J 
=P(J)«ST(J 
MUE 

1  =  1,4 
■H0 (I)»P(I 
=CO( I)*P(I 
=SO(I)#P(I 
)=HO( I)*PX 

)=co(  i)  »px 

)=Sl)(I)*PX 
)=HS(  D«OY 
) =CS ( T ) *0Y 
)  «SS  (  I  )  #qy 


)»XX-9. )*XM1/1A. 
1.125-P(2) 

• 

A. 

)#YY-«?.)*YM1/16. 

1.125-0(2) 

-0.125)«XP1 

(4) 

-1 .375)«xpl 

X  (?) 

-0.1?5)»YP1 

(4) 

-1.375)#YP1 
Y(?) 


)»KX  ♦  M  -KCOM 


tj) *HD( I) 
•J) ♦CD(T) 
•J) ♦SO(I) 
•I) *HS(I) 
tT)*CS(I) 
il) *SS( T) 


) ♦HP(K) 
) *CP(K) 
)+SP(K) 
( I) *HPX (K) 
(T) *CPX  (K) 
(I) +SPX (K) 
(I) +HPY (K) 
(I) *CPY(K) 
(I)*SPY(K) 


36 


?7  coimtinuf 
hso.o 

CKsU.O 
f^  =  0,0 

2H   MV=0.D 

hy=0.0 
CKX=0.0 

cky=o.o 

SKysO.O 
?*    NO  31  K=l,K4 

HrC(K)*HP(K) *H 
CK=C(K)»CP(K) *CK 
SK=C(K)*SP(K) *SK 

30  r<X  =  C(K)*HPX(K)*HX 
HY=C(K) »Hpy (K) *HY 
CKX  =  C(K)  <*CPX  (K)  *CKX 
C<Y=C (K)*C^Y (K) *CKY 
SKXsC  (K)#SPX  (K)  4-SKX 
SKY=C(K) *SPY (K) *SKY 

31  CONTINUE 

N4D=S0PTF (CK#CK*SK»SK) 
rKsCK/PAD 

sk=Sk/rao 

32  fJKX=CK*SKX«SK*CKX 

r)KY=CK»SKY-SK*CKY 
3^  KF.TUPN 

END 
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s 

c 
r 

c 

SPFKI) 

F 


UhRO 

OMMO 
'OMMO 
OMMO 
PA9 
ACl  = 
h  =  Fa 
*C1  = 
*C*  = 
VH  =  V 
4C1  = 
FsFfl 
AC1  = 
AC2s 
■Jh  =\/ 

AC1  = 
=  F  AC 
AC1  = 
Af?  = 
rt  =  B# 
=  (  VF 
C  =  A- 

a  =  ( n 

C  =  l)A 
X=Da 

YrUA 
V=D^ 
Y  =  U^ 
*=UC 
V=DC 

•in 


UTI 

N 

M/L 

M/L 

A!*F 

CI* 
2.1 
FAT 
H*F 

H  +  H 

ci* 

0.3 
FAC 
F»F 
H*3 
]  *■» 

3.H 

FAT 
FAT 
♦  \/H 

VH 

\/F  + 

-f)v 
»HX 

#Hy 
#MX 

*HY 

#hV 
*HV 
N 


Y.  A»R 
hX,HY 
A  X  ,  A  Y 
!vS  OF 
f]  088 
.  1  49h 
*94fr? 

0.0  93 
'? 

-»bAH6 

.34*3 

11219 

0.025 

? 

09733 

SM?43 

4  3  *  0  5 

0.0  hb 


MF. 

<  • 

?/ 

4/ 
7F 
n  . 
*2 

<*9 

ttC 

o 

• 

•  0 

4H 

1- 
AC 
7. 
^. 

»;? 

l- 

)  #0f5 


tCC«HtCK«SK»KX»KYtKTf FOXtFOYtKrON 

,HX,RY,CX,CY 

A     AP3 
«079 

94f>?5«*FXPF  (-0.093470S93044»H)  #0  .  0?534  1  0549*3 
5«/FACl 
470S93U44 

??74 

12l9fe96*F_XPF  ( -n. 025  324 954 73**H)»7. 963790 85?4 

696/FAC1 

324954736 

5?U 

60581 *F*PF  (-0,nHHo?^4S249ooH)  »0  ,40«48«7329l  F-4 

«l/FACl 

0?9452490 


l>  V  H )  #  0  •  5 

H 
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SUBROUTINE  VHERIV 

COMMON   X»Y«AtfltCCfH,CK»SK»KX»KY»KT»FOXtFOY»KrON 

C0MMON/L3/OKX»0KY 

COMMON/L4/AXfAY»BX»RY»CXfCY 

C0MM0N/L5/CAPViCAPVXtCAPVY»CAPVU.W 

REAL  MSEtMSFXtMSFYiNUM 

OELXsX-FOX 

l)FLY  =  Y-FOY 

MSF=(OEL  X*»2*OELY**2+973.75)/1043,63P743 

MSFx«OELx/52l .« 193715 

HSFY=0ELY/S21 .8193715 

TF(  LR  .EO.  0  )GO  TO  1 

COST  =  COSF(  W)*CK  ♦  SINF(ia/)*SK 

STNT  =  SINFM)«CK  -  COSF(W)»SK 

GO  TO  2 
1  COST=COSP*CK  ♦  SINP*SK 

STNT=SINP#CK  -  OnSP^SK 
?.    a2mC?=A»A-CC*CC 

POOTsSQRTF (B#»2* (COST»*2) +A2MC2* ( SINT«»2  )  ) 

NtlMsR#A2MC2/A 

DENOMbB*CC*COST/A*ROOT 

v=num/DEmom 

V<3RI0»V/B,  5660*16667 

CAPV=VGRIO*MSF 

AmCXsA*AX-CC*CX 

amcy»a*ax/cc*cy 
fx=(cc*bx*b*cx-b*cc*ax/a)/a 

FY= (CC*PY*fl*CY-R*CC*AY/A) /A 
BPCMA2sR»R  +  COCCA#A 
TERMX*(COST»FX*B*CC*SINT»DKX/A*(AMCX#SINT»SlNT*B»RX<»cOST»rOST    ♦ 

i  bpcma2*sint#cost«dkx)/root)»(a*vgrio/<b*a2mc2>    } 

tERMy3(C0ST*Fy*R*CC*SIK'T»DKy/A* (AMCY#SINT*SlMT*R*BY*rOST»COST 

1     ♦BPCMA2*SINT»cnST«DKY)/R00T)*(A#VORTD/(R»A2MC2>     ) 
VX=VGRID»(2.*AMCX/A2MC2    *BX/b-AX/a-TeRMX) 

v-ysVGRIH*  (2,*AMCy/A2MC2  +  By/B-Ay/A_TErMy) 

CAPVXaVX#MSF+VGRID*MSFX 

CftPVY=VY«MSF*VGRlO*MSFY 

VH=VGRID«(B*CC#SINT/A*8PCMA2»SINT»C0ST/R00T)/(B#CC#CoST/A    +ROOT) 

CAP'VU*MSF»VU 

UPTURN 

FND 
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SUPHOUTINiF    AMGLF 

COMMON)    X,Y»A,H»CC»H.CK»SKtKX,KYtKT«FnXtFOYtKCON 
COMMONl/l.ft/     XKiXLGtXLT 
ijFl  X  =  X-FOX 
:.)FLY  =  Y-FOY 

rnsXKs-.DFL.X*CK-nFLY*SK 
MMXK=OFl.X*SK-[)FLY»C* 
IF (COSXh )?»lf2 
1     kKsSTGMF(90«»MNxK) 

p  XKsATANF ( S T MXK/COSXK ) *5  7 . ?957795 1 
IF (COSXK) 3t**6 

3  1F(SIN>*)5»4,4 

4  XK=XK*1H0. 

00  TO  8 

b     xK=XK-]H(). 
h  [  F  ( X  K )  7  •  «  » ft 

h  xT=DFLX*.9ft4ftn775-DElY*. I  7364ft 18 
YTsUFI  x».l 73*481fl*DELY*.9848o775 
Pttp=SQPTF ( XT*XT     ♦     YTOYT) 

1  F  (  X  T  )  ]  0  ,  9  ,  1  0 

')     x|.  G  =  SI0N'F  (90.  «YT) 

So  to  J'* 
111  XLGsATANiF  (  YT/XT)  05  7.P9S7795J 

TF(XT)ll,14tl* 
11  I F  ( Y  T )  1 3 , 1 ?  »  1  ? 
\i>     XLGsXLG  ♦  1«0« 

.;?0  TO  14 
] 3  xLGsxLG  -  inn. 
j  <♦  XI  T  =  -ATAmF  (RAO/31  ,?f|C>)  #114. 59lSb9  *  ^O. 

PfTUPM 

pM|» 
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TDENT      DTGB 
PROGRAM    l.FMGTH 

BLOCKS 

wrOGram*    LOCAL 

ENTRY  POINTS 
nonouo  nTG8 

EXTERNAL  SYMHOLS 

TRFHTG 

FNTRY   DTGB 

«  FORTRAN  CALL  TO  THJS  PROGRAM-— CALL  DTGB(AtB,C) 

»  WHERE  A  ■  ADORERS  OF  OTG  WORD<LEFT  JUSTIFTEn)Bl 

»  Be  AOORE$S  TAy  B2 

*  C  *  ADDRESS  NEW  OTG  WORO(LEFT  JUSTIFlFD)B3 

»  RETURN  WITH  MEW  DTG  IN  C 

PTGn  f->Ss7    1 


RGT  ADJUST  DTG 


sx* 

B3 

SA* 

T 

SAi 

HI 

LXl 

4  8 

SA2 

B? 

PJ 

sXTRFOTG 

«;ai 

T 

«;h^ 

XI 

LX6 

12 

sa* 

B3 

JP 

OTGB 

DATA 

0 

END 

UNUSFO  STORAGE  2l  STATEMENTS         3   SYMBOLS 
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TDFMT  FIELOSX 
'POfiPnv    LENGTH 

BLOCKS 

3pnGRA^»         LnCAL 

"MTPY    POTmTS 

oonooo  fiflosx 
'xternal  symrols 

ryQF-i  T      i*AR64X  RDmST 

FN'TRY  FIELOSX 

»  REVISED  BY  S.W.  SF| FRIDGE  JUM  1968  MELLONICS 

» 

»           SELECT  FNWF  FIELD  FROM  MASTFP*  UNPACK  AND  FLOAT  DATA 

» 

»  FORTRAN  CALL  TO  THIS  PROGRAM—-  CALL  FIELDS  ( A  »PtCtD) 

»  WHEoF  A  b  ADDRESS  OF    39R9  WORD  P.IN  (CONTAINS  FLTD  DATA  ON  EXIT) 

»  R  s  ADDRESS  OF  DATF-TIMF  GROUP  WORD 

*  C  =  ADDRESS  OF  PaRa^ETfR  SELECTION  INngX <Q?1 t2»3) 
»  WHERE  0  3  CDtl  s  CHt2  s  U  CURR»3  B  v  CUPP 

»  D  b  ADDRESS  OF  PEAD-MASTER  ERROR  INDICATOR  WORD 

»  (If  fteLd  mot  fOUnd  or  Parity  error  (d>«1  on  return) 

fielosx  hss7  1 

»  Sawf  FORTRAN  CALL  PARAMETERS 

SA?  82                              OTG 

SX6  81 

&X7  X2 

c;A(6,  PAPAM 

SA7  TABLF+1 

SA3  R3                              STORE  SELFCTION  INDEX 

^X*  X3 

«;a*  counter 

sx7  r4                     not-found  indicator 

SA7  PARAM+1 
»          PFAD  ONE  FIELD  FROM  MASTER  TAPE 

BGM        SA2  COUNTER 

SAi  NAMTAR*x?                       BCD  PARAMFTER  NAME  OF  FIELD 

HXf,  XI 

^Afc  TABL.F 

SHI  A6 

PJ  bXRDMST 

•  UNPACK  FIELD 
<;ao  rO 
SHI  PAF 
SA?  PARAM 

SB?  X2                              LOC.  UNPACKED  FIELD 

SH3  7601P 

PJ  »XWAP64X 
•»          SHTFT  DATA  INTO  LOWER  4fl  BITS 

SB1  20 

SB?  3988 
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SB3 

BO 

SB7 

1 

SAl 

PARAM 

SHIFT 

SA2 

X1*B1 

AX2 

12 

RX6 

X2 

SA6 

Xl*B3 

*&\ 

Bl*B7 

SB3 

B3*B7 

GE 

B2iBl, SHIFT 

» 

CONVERT  FIELD  TO  FLOATING  POINT 

S81 

XI 

SB? 

XI 

SAl 

COUNTER 

SA? 

SCATAB*xl 

SB3 

X2 

SB4 

7601B 

RJ 

aXFXDFLT 

i 

JP 

FIELDSX 

TAP|  F 

DATA 

0 

DATA 

0 

vFn 

30/PAF,30/ERP 

DATA 

5LTAPE4 

DATA 

3LEEE 

NAMTAR 

DATA 

10HCD 

DATA 

10HCH 

DATA 

10HU   CURR 

DATA 

10HV   CURR 

SCAT  AH 

DATA 

41 

DATA 

41 

DATA 

36 

DATA 

36 

COUNTER 

DATA 

0 

PAPAM 

BSSZ 

2 

ERR 

SAl 

PARAM*1 

SB1 

1 

SX6 

Bl 

SAft 

XI 

JP 

FIELDSX 

PAF 

BSSZ 

EMO 

1100 

SCALE  FACtOP 


UNUSED  STORAGE 


85  STATEMENTS 


13   SYMBOLS 
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TOFNT   SEALAMO 
PPOGOAM  LENGTH 

MLOCKS 

PROGRAM*    LOCAL 

FMT^Y  POINTS 

on  no  (i  (i  sEalano 
c  x  rEpNAL  symbols 


I.ANDqffl 


SEA |  AMD 


St  A 

LAMP 

EWWO" 

OUT 


CMTRY   SEALAMD 

FOdTPAM  CALL  TO  THIS  PROGRAM CALL  SEALAND<If 

-WMfRE  I  =  ADDTSS  FO  RfaL  T  iNpEX  ON  FNWC  63Xf 

J  r  ApPRFSS  OF  REAL  J  INDEX 

i  AMn  =  AnnRcSS  r»r  I  amo/Sca  TNnTrATPR  t  n 


JtLANDtKEP) 
3  GRIo 


^SS7 
SA] 

SA? 

t  IX  1 

IX? 
SHl 
SH? 
I  xi 
1X2 

4X] 

RX? 

PX1 

SXfc 
SX7 
SA* 
SA7 
SHl 
SB? 
SH3 
JP 

MX* 
MX  7 
JP 
SX* 
MX  7 
JP 

SX* 
SX7 
SA] 
SA? 
SA6 
SA7 

JP 


J  r  APRRFbS  Oh  RFAL  J  INOEX 

LAND  =  AOOPESS  OF  LAND/SFA  INDICATOR  (0  EOR  SEAt  1  EOR 

KPR     s      APinDrC^     AC      rDDi-iD      TMnTrATAD      fOIIT     nr     aDIlMr^l 


keR  ■ 


AOpRESS  OF  ERROR  INDICATOR  (OUT 
(0  IF  ERROR  ARSENT.  1  IF  EPROR 


Of  bounds) 

IS  PRfSENT) 


LANf 


1 

PI 

P2 

XltBl 

X?»P? 

P1*3P 

R2*15 

XltBl 

x2,P2 

36 

xo*xi 

-x0#x2 

Xl*x2 

P3 

P4 

LSIM 

FRPIM 

SEA 

LANP 

ERROR 

sXLAMPSFA 

0 

n 

OUT 

1 

n 

OUT 

l 

1 

LSIN 

ERRIN 

XI 

X2 

SEA(_AND 

LSIN 
ERRIN 


I  S39 
J  $15 


I  BITS  24.47fJ  BITS  0-23t 
EACH  SCALFD  Sl5 

save  l/s  tnd.  address 
save  error  ind.  address 


MSS7 
WSSZ 
END 


DIMFMSloN  of  wavF  FIELD  aPRayS=(  22»  32t  12)     T'Uns  m 
ROUTE  OF  ^HIP  BEGINS    0  hours  AFTFP  TIME  SERIES  0«TGIN 
FROM  LONGITUDE  =   154.0  AN»n  LATITUDE  =    41.0 
TO  LONGITUDE  =  -123. n  ftNO  LATITUDE  =    3fl,0 


JMlMr  J? 


OF 


OAYC 

LOMGI 

1   A  T  J  - 

WAVE 

TRAV/EI 

-THOE 

TUOE 

HEIGH 

n  .on 

154.0 

4  1.D 

21.21 

•  5n 

lb6.6 

42.0 

19.31 

l.Oo 

159.9 

4  3 .  n 

16.76 

1 . 5  n 

163.6 

44,? 

15.63 

2.nn 

1*7.5 

45.  1 

20.  no 

2.5r 

171.0 

45.  H 

19.54 

3.00 

175.3 

46. b 

16.96 

3.5n 

It1?. 6 

47.  J 

11  .64 

4.00 

-1  /5.4 

47.5 

9.71 

4. bo 

-170.3 

47.7 

in. 36 

S.Oo 

-1*5.2 

47.7 

4.n9 

^.So 

-160.0 

<*7.4 

6.55 

6.  On 

-15S.0 

4  7.1' 

9.4ft 

^.Si 

-150.1 

46.3 

11.78 

7.00 

-145.5 

45.4 

17.28 

7. bo 

-141.2 

44.4 

13.41 

8.00 

-136.8 

43.2 

7.39 

B.5n 

-132.6 

41  .6 

4.ft6 

9.00 

-12«.b 

40.2 

1.91 

9. bo 

-124.6 

3*.S 

6.44 

9,7n 

-123.1 

37.8 

13.70 

(-fooestc  houte 

wave  direction 
from  modth 


31  Q 
304 
231 
252 
27o 
19<; 
266 
219 
109 
20M 
121 
232 
ISO 
175 
17* 
1*2 
153 
162 
29 
34* 
314 
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MINIMUM  TJMF  ROUTF 


DAY9 

OF    TRft\/F 

o, 

,0o 

i 

.4P 

.9a 

I 

.46 

1 

.9? 

2 

.40 

2 

,flp 

3, 

.3* 

3 

.84 

4, 

.3? 

4 

,8o 

5. 

.?n 

5, 

,7a 

ft, 

.24 

h 

.7? 

7, 

,?o 

7 

.ftp 

P 

,1ft 

8 

*ft4 

9, 

.1? 

9 

,ftn 

HE 

]        M       = 

LONG  I 

1  ATI 

-TH[)E 

TI'oE 

1^4,0 

41.0 

1S7.7 

40.4 

lhl.b 

4  0.3 

\b*,a 

4  0,5 

170.2 

40.9 

174.3 

41  .ft 

17*. 1 

42.7 

17H.0 

44.0 

173.8 

4S.1 

lft9.3 

45.9 

lft4.7 

4ft. 0 

lftO.O 

45. ft 

155.  ft 

44.9 

151.3 

^4.3 

147.3 

43.7 

14*. 4 

43.1 

H9.3 

42.5 

135.1 

4  1.7 

130.8 

40.7 

12ft. 8 

39.4 

123.0 

38.  1 

.00 


WAVE 

WA\/E    DIPEctTOn 

HEIGHT 

pROM    MQoTH 

21.21 

319 

14.80 

1 

ft. 92 

171 

10.10 

291 

14.70 

ISO 

IS.ftH 

281 

17.50 

2*0 

13.88 

275 

12.ft9 

2*2 

12. 4« 

275 

7.0ft 

U2 

6.«9 

3  34 

9.35 

147 

10.78 

157 

13.  «1 

lft4 

12.97 

171 

9.42 

lftft 

5.ft2 

2<*ft 

3.9ft 

IIS 

1.29 

1 

15.20 

P 

■•TIGHTEn    TTMF    ROUTF 


PET 


DAYc 

OF     TRAvF 

0, 

,00 

) 

.4* 

1 

.97 

1, 

,45 

1, 

,94 

2, 

.4? 

2, 

,9n 

3, 

,3o 

3, 

,A7 

4, 

.3* 

4( 

.84 

5, 

.3? 

5, 

.81 

ft, 

.29 

ft, 

.7o 

7, 

.2ft 

7, 

,74 

P, 

.23 

«, 

.71 

9, 

,2o 

9, 

.ftp 

LONG  I 
-TUDE 

154.0 
15  7.b 

lftl.ft 
1**,0 

170.2 
174.2 
177.9 

■17ft.  4 
174.3 
lft9.7 

■lft4.9 
lftO.l 
155.6 

•151.3 

•1*7,4 
143.5 
139.4 
13^.1 
130.9 

•12ft.  b 
122.9 


H      3 

I   ATI- 

T'nF. 

4  1.0 

39.9 
39. ft 

39.8 
4  0.4 
41.5 
42.9 
44.4 
45.9 
4ft.  8 
4f.O 
4ft. 4 
45.5 
44,7 
43. ft 
4  3.1 
42.4 
41.7 
4  0,7 
39.5 
3ft.  2 


.92 


WAVF 

WAV/f'   otpection 

HEIGHT 

pwOM     MOPTH 

21.21 

319 

13.39 

30 

5.18 

145 

9.10 

275 

13.90 

271 

15.39 

2P1 

17.81 

2P1 

13.48 

270 

11.71 

275 

11.14 

2*4 

5.44 

13ft 

8.41 

309 

9.7ft 

HP 

10.94 

158 

14.09 

Ift^ 

12.77 

171 

9.14 

lft7 

5.74 

1R7 

3.73 

109 

1.17 

359 

15.90 

1P1 
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GLOSSARY 
TEST 


IMIN   JMIN 

KX,  KY,  KT 

DTG(KT) 

U(200)  ,  DELU(200) 
GRADX,  GRADY,  GRADZ 
AZ,  AL ,  HT ,  WA 

WE(3989) ,  WD(63,63) 
XLAM,  XMU,  XSIGMA 

X,Y 

A,B,CC 

H,  CK,  SK 

FOX,  FOY 

XHT,  CSK,  SNK 
T 

TFS 

XMEAN,  YMEAN 
ROX,  ROY 


=  Minimum  indices  of  wave  heiqht 
data  subfield 

=  Dimensions  of  wave  data  sub- 
field 

=  Array  of  octal-date  time  groups 

=  Control  arrays 

=  Gradient  arrays 

=  Arrays  for  track  position 
tabulation 

=  Unpacking  area  of  core 

=  Components  of  adjoint  vectors 

12    3 

A  ,  A  ,  AJ 

=  Rectangular  coordinates  of  ship's 
position  relative  to  subfield 
origin 

=  Parameters  for  elliptical  polar 
velocity  diagram  Fig.  1 

=  Wave  height,  cosine  and  sine 
of  wave  direction  relative  to 
OXY  axes 

=  Floating  point  coordinates  of 
North  Pole  relative  to  subfield 
origin 

=  Wave  field  arrays 

=  Time  in  days  from  first  member 
of  time  series 

=  Time  in  days  from  beginning  of 
track 

=  Coordinates  of  midpoint  of  sub- 
field 

=  Semi-dimensions  of  a  rectangle 
in  the  subfield  within  which 
track  must  lie 
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LR 

C0SLG1,  SINLG1 
COSLTl,  SINLT1 
COSLG2,  SINLG2 
COSLT2,  SINLT2 
XI,  Yl,  X2,  Y2 

S12 
ARC 
XSTART,  YSTART 

XEND,  YEND 

COSP,  SINP 

V 
GAMA 


DELUMX 


N,  L 


=  Index  to  denote  a  geodesic  track 
(LR=0)  or  an  optimum  track 
(LR=1) 

=  Cosine  and  sine  of  longitude  of 
initial  point 

=  Cosine  and  sine  of  latitude  of 
initial  point 

=  Cosine  and  sine  of  longitude  of 
terminal  point 

=  Cosine  and  sine  of  latitude  of 
terminal  point 

=  Coordinates  of  initial  and  ter- 
minal points  of  track  relative 
to  North  Pole  in  grid  units 

=  Straight  line  distance  from  XI, 
Yl,  to  X2,  Y2 

=  Great  circle  distance  from  XI, 
Yl  to  X2,  Y2 

=  Coordinates  of  initial  point  in 
grid  units  relative  to  subfield 
origin 

=  Coordinates  of  terminal  point 
in  grid  units  relative  to  sub- 
field  origin 

=  Cosine  and  sine  of  ship's  head- 
ing at  any  point  along  great 
circle  track 

=  Control  program  to  define 
initial  control  u(t) 

=  Weighting  factor  used  in  cal- 
culating optimum  track,  GAMA=0 
minimum  time ,  GAMA  >  0  weighted 
time  track 

=  Maximum  allowable  change  in  the 
control  u(t)  at  anv  instant  of 
time 

=  Counters  to  stop  computation  if 
convergence  is  questionable 
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STEP 


XDOT,  YDOT,  ZDOT 


MSF 


All,  A12,  A22 


Zll,  Z12,  Z13,  Z22, 
Z23,  Z33 

DET 


CI,  C2 

El ,  E2 ,  DELTAU 

NN 


Bll,  B12,  B13,  B22, 
B23,  B33 


D 


Dl,  D2,  D3 


Bl,  B2,  B3 


RR 


=  Variable  time  steo  of  integra- 
tion 

=  Time  derivative  of  state 
variables 

=  Map  scale  factor  for  stereo- 
graphic  projection 

=  Elements  of  A  matrix  of  Sec- 
tion III 

=  Elements  of  Z  matrix  of  Sec- 
tion III 

=  Determinant  of  A  matrix  of 
Section  III 

=  Elements  of  C  matrix  of  Sec- 
tion III 

=  Elements  of  E  matrix  of  Sec- 
tion III 

=  Index  to  denote  a  minimum  time 
track  (NN=0)  or  a  weighted 
time  track  (NN=1) 

=  Elements  of  B  matrix  of  Sec- 
tion III 

=  Determinant  of  B  matrix 

=  Elements  of  D  matrix  of  Sec- 
tion III 

=  Elements  of  a  matrix  of  Sec- 
tion III 

=  Coefficient  of  DELU(T)  to  limit 
maximum  change  in  control  at 
any  time  to  DELUMX 


TERP 


C(4) 


TT 


=  Weighting  factors  for  interoola- 
tion  in  the  time  demension  be- 
tween K4  ordinates 

=  Index  to  determine  member  of 
field  time  series 

=  T  of  Eq.  (20)  of  reference  [5] 
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M,  N 
XX,  YY 
P(l)  to  P(4) 
Q(l)  to  Q(4) 
PX,  QY 


HX,  HY,  CKX,  CKY, 
SKX,  SKY 

DKX,  DKY 


VH,  VF 


AX,  BX,  CX,  AY,  BY,  CY 


MSFX,  MS FY 

COST,  SINT 

VGRID 

CAPV 

CAPVX,  CAPVY,  CAPVU 


=  Indices  to  pick  out  x,  v  grid 
point  data 

=  XY  of  equation  (18)  of 
reference  [4] 

=  Pi  to  P4  of  ecruation  (19)  of 
reference  [4] 

T 
=  Elements  of  P  (y)  matrix  of 

equation  (18)  of  reference  [4] 

=  Partial  derivatives  of  P(l)  to 
P(4)  and  0(1)  to  0(4) 

=  Partial  derivative  of  H,  cos  K, 
sin  K 

=  Partial  derivatives  of  wave 
direction  angle  K 


SHIP 


=  Ship's  speed  in  knots  in  head 
waves  and  following  waves 

=  Partial  derivatives  of  oaram- 
eters  of  elliptical  oolar 
velocitv  diagram 


VDERIV 


=  Partial  derivatives  of  map 
scale  factor 

=  Cosine  and  sine  of  the  angle 
0  of  Figure  1 

=  Ship's  speed  in  knots  relative 
to  the  stereographic  grid 

=  Ship's  speed  in  knots  relative 
to  the  earth's  surface 

=  Partial  derivatives  of  CAPV 
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